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Abstract 

The vibrated bed of powder, a vessel that mono disperse glass beads fill and a loud 
speaker shakes, is investigated numerically with the distinct element method, a kind 
of molecular dynamics. When the bed is heavily shaken, the displacement vectors 
of powder have the power spectrum with the dependence upon the wave number k 
as k~ 5 / 3 . The origin of this spectrum is suggested to be the balance between the 
injected and dissipative energy, analogous to the proposal by Kolmogorov to explain 
the fc -5 / 3 energy spectrum observed in the fluid turbulence. Furthermore, the same 
spectrum still appears even without flows of powder. Thus Kolmogorov's argument 
is more universal than believed before. 

Keywords: granular turbulence, Kolmogorov scaling, numerical simulation, vi- 
brated bed of powder. 
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1 Introduction 



For these a few decades, the non-linear physics, which has remained untouched over 
several years since there is not any convenient method to investigate it, has become 
one of the central topics in the modern physics. 

This is because the phenomenological approach provided by the statistical physics 
allows us to investigate them]!], 0. The basic philosophy underlying this approach 
is: 'Some phenomena are independent of the details of the specific systems; Hence, 
the phenomenological approach will be valid.' 

For example, the amplitude equation and the phase dynamics approach |l| can 
describe many phenomena ranging from the fluid dynamics to the chemical reaction. 
However, the validity of these approaches is restricted to the weak non-linear regime 
where only a few degrees of freedom survive. 

On the other hand, the systems where many degrees of freedom still survive 
and couples with each other — the subjects of the most interest among non-linear 
physicists — are attacked by the phenomenological numerical approaches, like the 
coupled map latticep|. Although they exhibit many interesting phenomena, they 
lack the direct relations to the real systems. At the moment, there is not any general 
tool to study the real non-linear systems with many degrees of freedom surviving. 

Instead of trying to find general tools, one can specify phenomena where many 
degrees of freedom survive, which appears in a real system with the easily treated 
model. It will become a start point to construct the general theory. 

An example among the phenomena that have the surviving many degrees of free- 
dom is the k~ 5 ^ 3 energy spectrum in the fluid turbulence, which is first proposed 
by Kolmogorov||. In his theory, he has predicted that the energy spectrum of the 
fluid turbulence, in the high enough wave number regime, has the dependence upon 
the wave number k as A; -5 / 3 , which was confirmed later experimentally and numer- 
ically. His theory is valid for almost all kinds of fluid turbulence:the experiments 
in the laboratories, the atmospheric flows — wind, and the flows in the tide. His 
success strongly supports the belief mentioned above:the details are not important. 
However, his theory is too general to make clear the mechanism from the dynam- 
ical points of view; it does not have any direct relationship to the Navier-Stokes 
equation:the basic equation of fluid motion. 

The attempts to solidify the foundation of his theory have not yet succeeded, 
since it is hard to obtain enough informations of the local structures of the fluid 
turbulence. Experimentally, the simultaneous measurements of the velocity over 
the whole region is impossible. Only variables observed experimentally is the time 
sequence of the velocity at a specific point, which is interpreted as the spatial struc- 
ture of the velocity by assuming unjustified Teylor's frozen turbulence hypothesis 0. 
Numerically, although the detail of the velocity field is easily obtained, the limited 
computational resources makes the production of fully developed fluid turbulence 
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difficult. 

Even phenomenological approaches are hardly applied to it. The amplitude 
equation or the phase dynamics approach cannot be applied to k~ 5 ^ 3 energy spec- 
trum where many degrees of freedom is essential. The phenomenological numerical 
approach has succeeded in reproducing them[||, |], but could not explain the mech- 
anism well due to the lack of the direct relationship to the experiment. Hence, for 
the fully understanding of the Kolmogorov's theory, it is hopeful to find the suitable 
substitution that obeys Kolmogorov's law and can be treated easily numerically and 
experimentally. 

In this paper, the k~ 5 ^ 3 energy spectrum is investigated in such a substitutiomthe 
granular flow. The granular systems, which many physicists are recently interested 
in|B|, |7|, ^ [|, are the systems whose number of degrees of freedom is limited enough 
to be integrated by the direct numerical simulations. For example, the number of 
degrees of freedom, even in the real experiments, is 10 6 , which should be compared 
with the Avogadro number 10 23 that represents those of typical materials like fluids, 
liquids, and solids. The numerical simulations of the system with one hundred 
degrees of freedom reproduce the non-trivial properties originating in the strong 
non-linearity as seen in the next section. Therefore, we can expect the experimental 
correspondence between the numerical model and the real experiments. 

The organization of this paper is as follows. In Sec.||, the results obtained nu- 
merically for the vibrated bed of powder is summarized. In Sec.^, to investigate the 
origin of the k~ 5 ^ 3 law, the new model that has the experimental correspondence is 
introduced. Discussions and conclusions will be presented in Sec||. 



2 Numerical Granular Turbulence in Vibrated Bed 
of Powder 



2.1 The experiment of vibrated bed of powder 

The experimental setup of vibrated bed of powder is as follows [|PI|, |l^|. The flat vessel 
with the horizontal dimension of about 10 cm is filled with granular matter:typically, 
mono disperse glass beads having the diameter of about less than 1 mm. The bed 
is shaken vertically with the acceleration comparable to the gravity. The speaker 
is usually used for this purpose, the frequency of the vibration is about from 10 to 
100 Hz, and the harmonic dependence upon time t is assumed;that is, the vertical 
displacement of the bed is as bcosuot. Although most experiments are performed 
in the three dimensional setups, the few did experiments in two dimensions. 

When the acceleration amplitude, T = bu^, exceeds some critical value T c slightly 
larger than gravity acceleration g, the bed exhibits several non-trivial dynamical in- 
stabilities:the heaping, the convection [|TT| , p~^ , IB] , |I4| , the capillarity ]Tj|, the surface 
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fluidization 
ing waves 



IB, [I7|], the Brazil nuts segregation |19| 20, |2T|. [2^], and the stand- 
and so on. 



23, 24, 25, 261 



2.2 Numerical Method — Distinct Element Method 

The distinct element method (DEM), or the particle element method, was named 
analogous to the finite element method. Unlike in the finite element method, each 
element in DEM can move freely, but like in the finite element method represents 
the macroscopic property of the material. The individual granular particle is treated 
as an element obeying the visco-elastic equation, 



N 
3=1 



k 




(1) 



where N is the total number of granular particles, Xj is the position vector of the 
ith particle, g is the gravity acceleration, and v is velocity. Each particle has a 
diameter of d. Because of a step function 6(x), particles interact only when they 
contact with each other, which represents the discrete nature of the granular matter, 
fco and r) are the elastic constant and the viscosity coefficient, respectively. This 
visco-elastic material constants give coefficient of restitution, e, and collision time 
t co i when two particles collide head-on; e = exp(— rjn/u) and t co i = tx/uj where u = 
(2k — rf) l l 2 . In this sense, the visco-elasticity is not the real material property, but 
the phenomenological one. The DEM can be also considered as a simple molecular 
dynamics simulation having the above interactions. 

Cundall and Struck[fZ7]], to investigate geological properties, has first proposed 
DEM which has come to have several versions later. However, in this paper, this 
simplified version, which can reproduce the convection |]T3"|, [T4|] — the most non-trivial 
phenomenon in the vibrated bed of powder, is used to study the vibrated bed of 
powder. 

The above equations are integrated by the Euler scheme where the time incre- 
ment At changes at each step so that the displacement during At of each particle 
does not exceeds some fixed value a (See appendix [A]). 

To simulate the vibrated bed of powder, the granular particles move within 
the two dimensional space — to save the cpu time-with the bottom oscillating as a 
function of time t, bcosuiot. In this representation the acceleration amplitude T, 
the strength of the external driven force, is hu^. When the bottom collides with 
a granular particle, it reflects the particle with the elastic constant k without the 
dissipation. 
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Figure 1: The process of fluidization. When T < r c , the bed consists of the solid 
region. As T exceeds T c , the fluid region appears, but the solid region still remains. 
For larger T, the bed is fully fluidized. 



Figure 2: Numerical setups used in this section. The bed is horizontally periodic 
with the period of L^, and the lid stays on the bed. 

2.3 Numerical results in the vibrated bed of powder 

First, the numerical behavior of the vibrated bed of powder obtained in the previous 



publications (T^, |28], |29[ is summarized. When T exceeds the gravity acceleration, 



the convection and the surface fluidization are observed ||13||. The fluidized region 
where the flow exists has the finite depth which increases as T increases and coexists 
with the solid region where the powder does not flow (See Fig.[I|). To understand 
the dynamics of the vibrated bed of powder, we consider the granular flow in detail. 

The granular flow is defined as the displacement of particle relative to each 
other. Thus, the motions like the translation and the expansion-not giving rise to 
the exchanges among particles-are not regarded as the flow in this subsection. 

The numerical setup to measure the flow is shown in Fig.0, while the two dimen- 
sional space is horizontally periodic and has the oscillating bottom. On the granular 
layer lies the lid which never tilts and does obey the equation of motion 

fh 
dt 2 



d 2 h N 

M— = -Y,0( yi -h)k o (h- yi ), (2) 



where M is the mass of the lid, h and yi are the vertical coordinate (height) of the 
lid and the ith particle, respectively. Each particle collides with the lid elastically — 
with k — and never pass through it. The lid suppresses the surface diffusion which 
causes the much larger displacement than the granular flow does. 

The lid is also necessary to know when to record the positions of particles to cal- 
culate the granular flow. When vibrating, the bed expands and contracts repeatedly 
and thus does not keep the constant volume. To remove the contribution of these 
motions to the displacements, which are much larger than those of the granular flow, 
we define the measuring time t n when the spacing Ah(t) (See Fig.|2|) between the lid 
and the bottom has some definite value, ho. It enables us to measure the displace- 
ment under the condition of the constant volume and to exclude the contribution of 
the volume changes to the displacement of each particle. Thus the definition of the 
displacement vector is as 

Ax| n) = Xj(t n+ i) - Xi(t n ). (3) 

The typical displacement vectors in the fully fluidized bed with T = 2.19 (u = 
2n/6,b = 2.0, N = 1024, L h = 128, e = 0.8,t coi = 0.1, d = 2.0, g = 1.0, M = 



5 



Figure 3: Typical displacement vectors of fully fluidized bed (T = 2.19) 

Figure 4: The dependence of the power spectrum S(K, Y) of the fully fluidized 
bed upon the wave number K and the height from the bottom Y (V = 2.19). (a) 
Y = 1 ~ 5,(b) Y = 6 ~ 10, (c) Y = 11 ~ 15, (d) Y = 16 ~ 20. The straight lines 
indicate the K~ 5 ^ 3 dependence. For large Y, i.e. near the surface, the spectrum 
becomes flat in the high wave number region. 

100, ho = 34.00) is shown in Fig.|| The spatial structure looks the velocity field in 
the turbulent fluid having the power spectrum with the k~ 5 ^ 3 dependence upon the 
wave number k. 

To compare it with the spatial structure of the velocity field in the fluid turbu- 
lence, Fourier power spectrum is calculated from the displacement vectors. First, I 
divide whole two dimensional space into d (horizontal) x y/3d/2 (vertical) cells. Each 
cell moves with the bottom and having numbers (X, Y) where X, Y are positive in- 
tegers. The cell (X, Y) covers the area d(X - 1/2) < x < d(X + 1/2), d(Y - 1/2) < 
y < d(Y + 1/2), where x and y are horizontal and vertical coordinates that have its 
origin on the bottom. 

I define displacement vector on each cell as: 

£ Ax?\ (4) 

ie(X,Y) 

where summation runs over only particles in the cell (X,Y). I calculate Fourier 
power spectrum for Ax^y on each layer 

S(K,Y) = (| Y.^x]y*M-3^XK/L) | 2 ) n , (5) 
x 

where K and L are integers, L is L^/d. Arr^y is the horizontal component of 

Ax^y. j is a pure imaginary number. The average (•••)« runs over 30 periods. 

Figure ^ shows the dependence of power spectrum upon the wave number K and 
the height from the bottom, Y . The power spectrum near the surface deviates from 
the straight line and has the flat spectrum — the white noise — in the higher wave 
number region, but for small Y, the power spectrum has the power dependence 
upon K. The slope of log- log plot is very close to -5/3 which coincides with that 
of the Kolmogorov's scaling theory. Therefore, we can conclude that the spatial 
structure in the vibrated bed is similar to that in the fluid turbulence. 

1 For the dense packing the particle forms the triangular lattice which has the width of Lh/d = 64 
and the hight of N/(Lhd) = 16. Thus the hight of the bed under the dense packing is ^ x d x 
{(N/L h d) - 1} - 26.0. 
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Figure 5: The dependence of the power spectrum S(K, Y) of the almost solidized 
bed upon the wave number K and the height from the bottom Y (T = 1.10). (a) 
Y = 1 ~ 5. (b) Y = 6 ~ 10. (c) y = 11 ~ 15. (d) Y = 16 ~ 18. For details, see 
Fig§ 

Figure 6: Typical displacement vectors of almost solidized bed. (r = 1.10) 

The origin of this turbulent motion in the fluidized region should exist in the 
motion of solid region, because as shown in Fig.|l| the solid region becomes fluidized 
gradually as T increases. To compare the motion in the solid region with that in 
the fluidized region, the power spectrum averaged over 64 periods in the solid region 
when T = 1.10(6 = 1.0, ho = 30.0, and the remaining parameters are identical to 
the above) is shown in Fig.[5[ As seen in the typical displacement vector (Fig||), 
almost all region is solid — no exchange between particles, but the global structure of 
the power spectrum is qualitatively similar;the power spectrum in the lower layers 
has the power dependence upon K with the exponent —5/3, and the white noise 
appears in the higher wave number region for the upper granular layers. 

In the following, I regard this k~ 5 ^ 3 spectrum without flow as the origin of that 
seen in the fluidized region and make clear how the k~ 5 ^ 3 spectrum appears without 
the flow. 

3 The /c~ 5 / 3 power spectrum without flow 

To understand the mechanism of the k~ 5 ^ 3 power spectrum without flow, a model 
in which the powder does not flow is proposed. Each element, the granular parti- 
cle, forms a triangular lattice?] as shown in Fig.|2] and interact, visco-elastically as 
described in Eq. ([!]) but without gravity g, with the six neighbors. The bottom oscil- 
lates as before, but the lid does not obey the equation of motion, Eq.(||). Instead, the 
lid oscillate as a function of time, b cosu Q t. That is, when the vertical coordinate 
of the bottom oscillate as b cos u t, then the height of the lid h = h + b Q cos u t. 
Thus, ho becomes the average distance between the bottom and the lid. At the 
initial stage, the granular particle lies on the triangular lattice between the bottom 
and the lid and moves obeying the Eq.(fJ) with g = 0. 

This numerical setup corresponds to the horizontally vibrated bed in the real 
experiments. One prepares the flat cell filled with the glass beads and makes two 
side walls oscillate harmonically. Hence, the following results are expected to be 
reproduced in the experiments. 

2 Hence no granular flow — no exchange between particles — is allowed at all 
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Figure 7: The time development of the energy spectrum, E(K).(a) The initial 
state(t = 1.5) and the early stage. (b) The medium stage where the high wave num- 
ber components start to have flat part. (c)The late stage. The energy spectrum 
shows the system reach the thermal equilibrium state. 

First, both the lid and the bottom oscillate with the same harmonic form, i.e. 
& = &o = l-0,g = 0.0. The ho is taken to be equivalent to the width of the triangular 
lattice having N/ L rows, \f?>/ 2 x \{N/ L) — 1] , and thus the granular particle is packed 
without space between each other. The remaining parameters are the same as used 
in the previous section. 

Figure [7] shows the time development of the energy spectrum 

L 

E{K,t) = (| ^Av x {t;iJ)exp(-j2mK/L) \ 2 )j, in mai, (6) 
i=i 

where v x (t; i,j) is the x component — the direction perpendicular to the direction of 
the vibration — of the velocity at the site at time t. The average is taken over 
the layers — y direction — and the ten initial realization of the velocity with the white 
noise energy spectrum. One should notice that it is the true energy spectrum — not 
the power spectrum of the displacement vector as in the previous section — thus 
we can compare it with the Kolmogorov's argument. At the initial time t = 1.5, 
each granular particle has the small random velocities and thus has the flat energy 
spectrum — the white noise. As time proceeds, the energy spectrum starts to incline 
and comes to have power dependence upon the wave number k. Furthermore, its 
slope is very close to the -5/3. However, for the medium stage, the flat spectrum 
appears in the higher wave number region and reaches the lower wave number region. 
Finally in the later stage, the total energy spectrum has become flat. 

This process can be interpreted as follows. First, the energy starts to dissipate in 
the high wave number region where the energy dissipation rate takes the large value. 
It enables the system to reach the stationary state described by the Kolmogorov 
theory. However, the further injection of the energy dominates the dissipation and 
the system reaches the thermal equilibrium state:each mode has the same amount of 
the kinetic energy. This is also seen in the numerical simulation of the Navier-Stokes 
turbulence when the dissipation is not large enough. Therefore, the dynamics in the 
vibrated bed of powder is essentially equal to that of the fluid turbulence. 

However, the k~ 5 ^ 3 spectrum appears only temporally, contrast to the permanent 
appearance of it in the previous section. The dissipation in the vertically vibrated 
bed may be larger because the bed is compressed between the lid and the bottom. 
To take into account this effect, the value of b is taken to be while &o remains 
unchanged so that the volume of the lattice changes. Starting from the initial white 
noise, the power spectrum becomes k~ 5 ^ 3 one and returns to the flat one. However, 
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Figure 8: The dependence of the energy spectrum E(K) upon the fraction of the 
cycle, toot = nn/6.(n = 1, 2, . . . , 6). The k~ 5 ^ 3 dependence is shown for the compar- 
ison. 

Figure 9: the total energy, the energy input rate, and the energy dissipation rate 
as a function of time t (See appendix |B[). 

this time the power spectrum recovers k~ 5 ^ 3 dependence from the white noise and 
this process is repeated once a period of the oscillation(Fig.|^). Here the energy 
spectrum is averaged over the ten periods to show the dependence upon the fraction 
of cycle. Thus k~ 5 ^ 3 power spectrum can appear repeatedly if the effect of the 
compression is considered. 

4 Discussion and Summary 

In the previous section, the powder bed has the energy spectrum with k~ 5 ^ 3 . The 
results also suggest that the dissipation plays the important role. The behavior of the 
energy spectrum is coincident with that in the Navier-Stokes numerical simulation 
when the dissipation is not large enough. However, there is not the flow at all, thus 
one cannot expect the powder obeys the Navier-Stokes equation. So, what is the 
relation to the Kolmogorov theory proposed for the fluid flow? 

Here the Kolmogorov's argument should be reconsidered. The essential assump- 
tion in the Kolmogorov's theory is that the balance between the energy input and 
the energy dissipation governs the dynamics of system. This means, the amount of 
total energy is as large as the product of the energy dissipation rate — or the energy 
input rate — and the characteristic time length. Figure [5] shows the total energy, the 
energy input rate, and the energy dissipation rate in the fluidized region investigated 
in Sec.0 as a function of time t. The total energy consists of both the elastic energy — 
among particles, the bottom, and the lid — and the kinetic and gravitational energy 
of particles and the lid. The dissipation energy comes from the viscosity between 
particles, and the work done by the vibrating bottom is considered as the energy 
input. The product between the energy input rate (~ a few thousands) and the 
period of the vibration (=6) is about 10 4 . On the other hand, in Fig|| the actual 
total energy is about 10 4 , when we eliminate the gravitational potential of the grand 
state(~ 2 x 10 4 ). Thus, the total energy included in the bed is almost equal to the 
energy injected each period. This means, the balance between the energy input 
and the energy dissipation dominates the dynamics of the system. So far, the basic 
assumption necessary for the Kolmogorov's theory is satisfied. 

The energy transport between the modes with the different wave number is also 
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Figure 10: The schematics of the energy flow in the vibrated bed. 



assumed to be caused by the non-linearity, which the step function 6(x) in Eq.(JT]) 
has. Actually speaking, Kolmogorov has never used the existence of the flow to 
derive his scaling form, since it is apparent;the turbulence must flow. Here, the 
steady state which he proposed in his scaling theory can appear without the flow 
because the vibrated bed of powder can satisfy the basic requirements without the 
flow. 

Another difference from the fluid is the dimensionality. The A; -5 / 3 spectrum can 
appear only in the three dimension, not in the two dimension where it appears in 
the powder. However, again, Kolgomorov has not used the dimensionality explicitly. 
The main difference between two and three dimension is whether there is another 
conserved variable, the enstrophy^]. The enstrophy is impossible to define in the 
vibrated bed of powder because in the discrete material like the powder the spatial 
derivative, which is necessary to calculate the enstrophy, does not exists due to the 
lack of the continuity. Hence the enstrophy, which does not exist, cannot prevent 
the two dimensional system of the powder from obeying Kolmogorov's argument. 

Assuming the /c _5//3 power spectrum comes from the Kolmogorov theory, we can 
explain the power spectrum seen in Sec. |[ For lower layers of the bed, the compres- 
sion is strong enough, thus the clear k~ 5 ^ 3 spectrum can be observed. However, the 
upper layers are not compressed enough to dissipate injected energy, which makes 
the energy stay in the high wave number region and causes the flat spectrum. Figure 
|10| shows the schematics of the energy flow in the vibrated bed considering the above 
discussions. 

In summary, the vibrated bed of powder has the k~ 5 ^ 3 power spectrum of the 
displacement vector. It comes from the local motions in the solid region which the 
Kolmogorov theory can explain. This means, I found that the steady state proposed 
by the Kolgomorov can appear even without flow;thus his theory is more universal 
than believed before. The real experiments should confirm these results. 

This finding will make the effort to understand the k~ 5 ^ 3 mechanism easier, since 
the powder system does not flow — thus we can use lattice representation, and the 
dimensionality is two — requiring smaller computational resources than the three 
dimension does. The efforts also may give rise to understand the general feature 
of the strongly non-linear system where many degrees of freedom survives and may 
give us the general method to study such systems. 



3 The enstrophy is defined as (w(x, t) 2 }/2, where u> is the vorticity, rotv. 
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A Numerical method to integrate the equation of 
motion 

The method employed to integrate the equation of the motion Eq.(|l|) is Euler scheme, 

x 4 (t + At) = Xi(t) + v 4 (t)At + ^a,(t)(At) 2 (7) 
v 4 (t + At) = Vi (t) + £k(t)At, (8) 

where aj is the acceleration of ith particle. Eq.(^) enables us to relate the time 
increment At to Axi(t) =| Xi(t + At) — Xi(t) | and Ay^t) =\ yi(t + At) — y^t) |, 
where Xi and yi are the two components of x, respectively. Then 



Jv±i + 2a xi Axi — v xi 
A«„ = l-H 1 (9 ) 



At tl - ,* + (10) 

(11) 

where v X i and v y i are the x and components of Vj, and a x i and are the x and 
y components of a«. Using these equations, the representation of At with Axj, At 
takes the value of min Xi At Xi with the fixed Axi and Ayi. In the present simulation, 
Axi = Ayi = a = 0.005d. 



B The calculation of the energy and dissipation 

In this appendix, the definition of several quantities appearing in Fig.|9| is explained. 

The total energy consists of the kinetic energy, the elastic energy, and the grav- 
itational potential energy. The kinetic energy E kinetic is defined as, 

1 N 1 

E kinetic = -Mv 2 d + £ -v, 2 , (12) 

A i=i A 
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where Vud is the velocity of the lid. The gravitational potential energy U g can be 
expressed easily, too, 

N 

u g = Mgh + J2gyi- (is) 

i=i 

Here one should remember h is the height — the vertical coordinate — of the lid. The 
elastic energy U et consists of those among the particles, the bottom plate, and the 
lid, 



U ei = \k l^20(d- | Xj - x,- \)(d- | Xj - x 



2 

N 



■3 I)' 



+ ^6»(6cosw t - Ui)(b cos u t - ytf 



i=i 



iv i 

+ ^^-^(y,-^ , (14) 
i=l J 

Thus the total energy E tot can be expressed as E tot = E kinetic + U g + U e i. 

The energy input rate E input is the amount of work done by the bottom per unit 
time, 

1 N 

E in P ut = -r- J2i k o { e ( bcos ^ot-yi{t)){b cos uj t-yi(t))} 

x5{coscu (t + At) - cosujt}}. (15) 

Typically, At is taken to be 1 x 1(T 5 . The energy dissipation Edi S due to the viscosity 
between particles % and j is rj f | (vj — Vj) ■ rf(xj — x,) | when two particles contact 
with each other. Using d(xj — ~x.j)/dt = (vj — Vj), we get = ?7 /(vj — \j) 2 dt. In 
the present discretization substituting Vj(t + At) = Vj(t) + aj(t)At, 

^ = xt ? 1 x * " Xj l} { Av ?> At + AVij ' a ^ (At)2 + l a U At ) 3 } > ( 16 ) 

where Av tJ = Vj — Vj and Aa^ = a^ — a.j. 

These above variables have the relation, E tot (t + At) = E tot (t) + {E input (t) — 
E dts (t)}At + 0(At 2 ). 
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